! VSemiletov, 07.2011
! To set cylinder

module Acoustics
    ! SurfaceData

    logical :: AcousticSurf=.true.  ! on/off
    logical :: LeftDiscActive=.true., RightDiscActive=.true. ! on/off
    integer :: nPoints = 0          ! number of countour points
    integer :: nMaxPoints = 2000
	real(8) :: Contour3D(3,2000)    ! 3D coord of convex contour
    real(8) :: Center(3)            ! coord Center of convex contour // left disc's center
    real(8) :: x(3),y(3)            ! 2D basis in plane of the convex contour
    real(8) :: Contour2D(2,2000)    ! 2D coord of the convex contour in 2D basis
! input here! NOTE, Semiletov
    ! radius-vector from th left disc's center to the right disc's center
    real(8) :: CenterDir(3)=(/0.,1.0,0./)
    real(8) :: Length = 0.2 ! cylinder length
    logical :: Sphere = .false.
    real(8) :: Radius = 1d0
    ! factor of contour decreasing from the right side
    real(8) :: factor=1.

    type tAcousticSurface
        character(280) :: name
        integer :: nameLength
        character(280) :: fileContourIn    ! the name of the file for contour reading
        character(280) :: fileContourOut    ! the name of the file for contour reading
        logical :: DefCellsFromContour = .true. ! .true. = calc. Cells from contour, else read Cells from file
        integer :: nCellsOut = 0 ! number of cells for output data
        integer :: nCellsOutGlobal = 0 ! number of cells for output data in all processors
        integer,allocatable :: CellsOut(:) ! number of cell for output
        real(8),allocatable :: Area(:)
        real(8),allocatable :: xOut(:,:)
        real(8),allocatable :: normal(:,:)
        
        logical :: CalcAcoustics = .true.
        character(50) :: AcousticsOut
        integer :: LengthAcousticsOut

        ! Contour Definition:
        logical :: LeftDiscActive=.false., RightDiscActive=.false. ! on/off
        real(8) :: CenterDir(3)=(/0.,1.0,0./)
        real(8) :: Length = 0.2 ! cylinder length
        real(8) :: factor = 1.0d0

        !Sphere
        logical :: Sphere = .false.
        real(8) :: Center(3) = 0d0
	    real(8) :: Radius = 1d0

        real(8) :: timeStart=0.0d0, timeEnd=1.0d10  ! time interval for output
        real(8) :: dTimeOut = 1d-11 ! how often to output data
        real(8) :: nextTimeOut = 0d0 ! when acoustic output needs to be done
        integer :: ndTimeOut = 0 ! number of output times
        
       	real(8),allocatable :: recv(:)               ! for sending Volume between Points from diff procs
        integer,allocatable :: recvCounts(:)
        integer,allocatable :: rDispls(:)
        
        logical :: MeshDataOut = .true. ! output like hsp
        
!        logical :: DefCellsFromboundary = .false.
        integer :: nBoundary = 0
    end type tAcousticSurface

    integer :: nAcousticSurfaces=0 ! number of Acoustic Surfaces    
    type(tAcousticSurface),allocatable :: AcousticSurface(:)

    ! ListenerData
    type tListener
        character(280) :: name
        integer :: nameLength
        real(8) :: Location(3)=0.0
    end type tListener
    
    integer :: nListeners=0 ! number of Listeners
    type(tListener),allocatable :: Listener(:)

    ! Noise Signal
    type tNoiseSignal
        integer :: iSurface=0  ! number of Surface
        integer :: iListener=0 ! number of Listener
        integer :: nameFileLength=0
        character(280) :: outputFile    ! the name of the file for output noise signal
        
        integer :: nValues=0 ! number of time's points for output
        integer :: nValuesGlobal ! max number of time's points for output from all processors
        real(8) :: Pobs(300000)=0 ! pertubation of pressure
    end type tNoiseSignal

    integer :: nNoiseSignals=0
    type(tNoiseSignal),allocatable :: NoiseSignal(:)

contains

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! HeavysideFunc(r) = 1,   if r is outside of the surface
! HeavysideFunc(r) = -1,  if r is inside  of the surface
! HeavysideFunc(r) = 0,   if r is on the side where disc is not active
integer function HeavysideFunc(r)
!use Acoustics, only: AcousticSurf, LeftDiscActive, RightDiscActive, &
!                                nPoints, Contour3D, Center, x,y, Contour2D, &
!                                CenterDir, Length, factor
implicit none
    real(8),intent(in) :: r(3)
    
    integer :: iPoint
    real(8) :: ScalarProd, localFactor, vec3D(3)
    real(8) :: vec2D(2), SurfVec(2), SurfNormal(2)
    real(8) :: ScalarProdCenter, rMod
    
    HeavysideFunc = 0
    
    vec3D(1:3) = r(1:3) - Center(1:3)

! spherical surface:
    if(Sphere)then
        rMod = sqrt(sum(vec3D(1:3)*vec3D(1:3)))
        if(rMod.lt.Radius)then
            HeavysideFunc = -1
            return
        else
            HeavysideFunc = 1
            return
        end if
    end if

! from the left side of cylinder
    ScalarProd = sum(vec3D(1:3)*CenterDir(1:3))
    if(ScalarProd.lt.0)then
        if(LeftDiscActive)then
            HeavysideFunc = 1
        end if
        return
    end if
    
! from the right side of cylinder
    if(ScalarProd.gt.Length)then
        if(RightDiscActive)then
            HeavysideFunc = 1
        end if
        return
    end if
    
    localFactor = 1 + (factor-1.) * ScalarProd / Length
    vec2D(1) = sum(x(1:3)*vec3D(1:3)) / localFactor
    vec2D(2) = sum(y(1:3)*vec3D(1:3)) / localFactor
    
    do iPoint = 1,nPoints-1
        SurfVec(1:2) = Contour2D(1:2,iPoint) - Contour2D(1:2,iPoint+1)
        SurfNormal(1) = - SurfVec(2)
        SurfNormal(2) = SurfVec(1)
        
        ScalarProdCenter = sum(SurfNormal(1:2)*Contour2D(1:2,iPoint))
        ScalarProd = sum(SurfNormal(1:2)*(vec2D(1:2)-Contour2D(1:2,iPoint)))
        
        if(sign(1.0d0,ScalarProdCenter).eq.sign(1.0d0,ScalarProd))then
        ! 1. -> 1.0d0     SK: needs to be real for some compilers 
            HeavysideFunc = 1
            return
        end if
    end do

    HeavysideFunc = -1

end function HeavysideFunc
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

subroutine DefNormalDirection(xCenter,meshSf)
implicit none
    real(8),intent(in) :: xCenter(3)
    real(8),intent(inout) :: meshSf(3)

    real(8) :: xCenter3D(3)
    real(8) :: xCenter2D(2), meshSf2D(2), r(2)
    real(8) :: scalarProd

    xCenter3D(1:3) = xCenter(1:3) - (Center(1:3) + 0.5d0 * Length * CenterDir(1:3))

    scalarProd = sum(xCenter3D(1:3) * meshSf(1:3))

    if(scalarProd.lt.0)then
        meshSf(1:3) = -meshSf(1:3)
    end if

end subroutine
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
end module Acoustics
